
load("RData.sims.B.nobs1000")

#source("pscore_estimation_core1.R")
library(Matching)

nboots <- 0
print.level <- 0

#load(file="RData.simdata.B.obs1000.reps1000")
#dta <- simdata.B.obs1000.reps1000
#rm(simdata.B.obs1000.reps1000)
#gc()

sims.run <- 1000

logit.bm1 <- list()
boost.bm1 <- list()
forest.bm1 <- list()
nnet.bm1 <- list()
genmatch.bm1 <- list()

genmatch.weights <- read.csv(file="simsB1.nobs1000.gm1.results", header=FALSE)

for (s in 1:sims.run)
  {
    cat("\n\ns:",s,"\n")

    sdta <- as.data.frame(dta[,s])

    #true model
#    lm1 <- lm(y.a~ z.a+ w1 + w2 + w3 + w4 + w8 + w9 + w10, data=sdta)
#    true.est[s] <- lm1$coef[2]
#    cat("truth:",mean(true.est[1:s]),"\n")

    #logit
    logit.m1 <- Match(Y=sdta$y.a, Tr=sdta$z.a, X=logit.pscore[[s]], estimand="ATT", M=1)
    summary(logit.m1)

    logit.bm1[[s]] <- MatchBalance(z.a ~ w1+w2+w3+w4, match.out=logit.m1, data=sdta, nboots=nboots, print.level=print.level)
#    logit.bm2 <- MatchBalance(z.a ~ w5+w6+w7, match.out=logit.m1, data=sdta, nboots=nboots, print.level=print.level)
#    logit.bm3 <- MatchBalance(z.a ~ w8+w9+w10, match.out=logit.m1, data=sdta, nboots=nboots, print.level=print.level)    

    #boosting!
    boost.m1 <- Match(Y=sdta$y.a, Tr=sdta$z.a, X=boosting.pscore[[s]], estimand="ATT", M=1)
    summary(boost.m1)

    boost.bm1[[s]] <- MatchBalance(z.a ~ w1+w2+w3+w4, match.out=boost.m1, data=sdta, nboots=nboots, print.level=print.level)
#    boost.bm2 <- MatchBalance(z.a ~ w5+w6+w7, match.out=boost.m1, data=sdta, nboots=nboots, print.level=print.level)
#    boost.bm3 <- MatchBalance(z.a ~ w8+w9+w10, match.out=boost.m1, data=sdta, nboots=nboots, print.level=print.level)        


    #Random Forest
    forest.m1 <- Match(Y=sdta$y.a, Tr=sdta$z.a, X=forest.pscore[[s]], estimand="ATT", M=1)
    summary(forest.m1)

    forest.bm1[[s]] <- MatchBalance(z.a ~ w1+w2+w3+w4, match.out=forest.m1, data=sdta, nboots=nboots, print.level=print.level)
#    forest.bm2 <- MatchBalance(z.a ~ w5+w6+w7, match.out=forest.m1, data=sdta, nboots=nboots, print.level=print.level)
#    forest.bm3 <- MatchBalance(z.a ~ w8+w9+w10, match.out=forest.m1, data=sdta, nboots=nboots, print.level=print.level)            


    #nnet
    nnet.m1 <- Match(Y=sdta$y.a, Tr=sdta$z.a, X=nnet.pscore[[s]], estimand="ATT", M=1)
    summary(nnet.m1)

    nnet.bm1[[s]] <- MatchBalance(z.a ~ w1+w2+w3+w4, match.out=nnet.m1, data=sdta, nboots=nboots, print.level=print.level)
#    nnet.bm2 <- MatchBalance(z.a ~ w5+w6+w7, match.out=nnet.m1, data=sdta, nboots=nboots, print.level=print.level)
#    nnet.bm3 <- MatchBalance(z.a ~ w8+w9+w10, match.out=nnet.m1, data=sdta, nboots=nboots, print.level=print.level)


    #GenMatch
#    s <- 645
#    sdta <- as.data.frame(dta[,s])    
    Xm <- cbind(logit.pscore[[s]], sdta$w1, sdta$w2, sdta$w3, sdta$w4, sdta$w5, sdta$w6, sdta$w7, sdta$w8, sdta$w9, sdta$w10)    
    genmatch.m1 <- Match(Y=sdta$y.a, Tr=sdta$z.a, X=Xm, estimand="ATT", M=1, Weight.matrix=diag(as.real(genmatch.weights[s,])))
    summary(genmatch.m1)
    print(genmatch.m1$est)

    genmatch.bm1[[s]] <- MatchBalance(z.a ~ w1+w2+w3+w4, match.out=genmatch.m1, data=sdta, nboots=nboots, print.level=print.level)
#    genmatch.bm2 <- MatchBalance(z.a ~ w5+w6+w7, match.out=genmatch.m1, data=sdta, nboots=nboots, print.level=print.level)
#    genmatch.bm3 <- MatchBalance(z.a ~ w8+w9+w10, match.out=genmatch.m1, data=sdta, nboots=nboots, print.level=print.level)    
  }#end of sims loop

logit.pvals <- matrix(0, nrow=sims.run, ncol=1)
boost.pvals <- matrix(0, nrow=sims.run, ncol=1)
forest.pvals <- matrix(0, nrow=sims.run, ncol=1)
nnet.pvals <- matrix(0, nrow=sims.run, ncol=1)
genmatch.pvals <- matrix(0, nrow=sims.run, ncol=1)

for (ii in 1:sims.run)
  {

    logit.pvals[ii] <- logit.bm1[[ii]]$AMsmallest.p.value

    boost.pvals[ii] <- boost.bm1[[ii]]$AMsmallest.p.value

    forest.pvals[ii] <- forest.bm1[[ii]]$AMsmallest.p.value

    nnet.pvals[ii] <- nnet.bm1[[ii]]$AMsmallest.p.value

    genmatch.pvals[ii] <- genmatch.bm1[[ii]]$AMsmallest.p.value


  }

save.image(file="balance3.simsB1.nobs1000.Rdata")

box.data <- as.real(c(logit.pvals, boost.pvals, forest.pvals, nnet.pvals, genmatch.pvals))
type <- c(rep(1,sims.run), rep(2,sims.run), rep(3,sims.run), rep(4,sims.run),rep(5,sims.run))

box.data <- as.data.frame(cbind(as.real(box.data), as.factor(type)))

pdf(file="balance3.simsB1.nobs1000.pdf")
boxplot(V1~V2, data=box.data, names=c("logit","BOOST","RFRST","NN","GenMatch"), ylab="p-values")
dev.off()

#quartz()
#boxplot(V1~V2, data=box.data, names=c("logit","BOOST","RFRST","NN","GenMatch"), ylab="p-values")




